library(survey)
library(foreign)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(dplyr)

## see PCM_cleaning.R for data cleaning
load("~/Dropbox/Public_Conf_Mil/Data_and_code/wave1_data/clean_dataw1.RData")

## Create weighted survey design object
milconf_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

#### MILITARY ATTRIBUTES BY PARTY ####

## Most members of the military are educated - Q27A
svyby(~ Q27A_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27A_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are selfless - Q27B
svyby(~ Q27B_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27B_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are homophobic - Q27C
svyby(~ Q27C_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27C_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are racist - Q27D
svyby(~ Q27D_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27D_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are heroes - Q27E
svyby(~ Q27E_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27E_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military have PTSD - Q27F
svyby(~ Q27F_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27F_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are truthful - Q27G
svyby(~ Q27G_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27G_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are dangerous - Q27H
svyby(~ Q27H_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27H_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are sexist - Q27I
svyby(~ Q27I_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27I_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

## Most members of the military are intelligent - Q27J
svyby(~ Q27J_fac, ~ party, milconf_design, svymean, na.rm = TRUE)
svyby(~ Q27J_fac, ~ party+DOV_ASSIGNMENT_A, milconf_design, svymean, na.rm = TRUE)

#### FIGURE 2.8 - DEMOGRAPHIC FEATURES ####

dm_prty <- as.matrix(svyby(~ Q8_mr, ~ party, milconf_design, svymean, na.rm = TRUE))
dm_prty <- cbind(c(rep(15,3)), dm_prty)
dm_prty[,2] <- c("Democrat", "Independent", "Republican")

dm_ideo <- as.matrix(svyby(~ Q8_mr, ~ ideo3, milconf_design, svymean, na.rm = TRUE))
dm_ideo <- cbind(c(rep(14,3)), dm_ideo)
dm_ideo[,2] <- c("Liberal", "Moderate", "Conservative")

dm_gend <- as.matrix(svyby(~ Q8_mr, ~ male, milconf_design, svymean, na.rm = TRUE))
dm_gend <- cbind(c(rep(13,2)), dm_gend)
dm_gend[,2] <- c("Female", "Male")

dm_race <- as.matrix(svyby(~ Q8_mr, ~ RACETHNICITY, milconf_design, svymean, na.rm = TRUE))
dm_race <- cbind(c(rep(8,4)), dm_race[c(1, 2, 4, 6),])
dm_race[,2] <- c("White", "Black", "Hispanic", "Asian")

dm_relg <- as.matrix(svyby(~ Q8_mr, ~ religion, milconf_design, svymean, na.rm = TRUE))
dm_relg <- cbind(c(rep(9,4)), dm_relg[2:5,])
dm_relg[,2] <- c("Christian", "Catholic", "Other", "None")

dm_mvet <- as.matrix(svyby(~ Q8_mr, ~ vetstatus, milconf_design, svymean, na.rm = TRUE))
dm_mvet <- cbind(c(rep(12,3)), dm_mvet)
dm_mvet[,2] <- c("Active Duty", "Veteran", "Civilian")

dm_soco <- as.matrix(svyby(~ Q8_mr, ~ social, milconf_design, svymean, na.rm = TRUE))
dm_soco <- cbind(c(rep(10,2)), dm_soco)
dm_soco[,2] <- c("No", "Yes")

dm_fami <- as.matrix(svyby(~ Q8_mr, ~ family, milconf_design, svymean, na.rm = TRUE))
dm_fami <- cbind(c(rep(11,2)), dm_fami)
dm_fami[,2] <- c("No", "Yes")

dm_genr <- as.matrix(svyby(~ Q8_mr, ~ generation, milconf_design, svymean, na.rm = TRUE))
dm_genr <- cbind(c(rep(6,5)), dm_genr)
dm_genr[,2] <- c("Silent", "Boomer", "Gen X", "Millennial", "Gen Z")

dm_regn <- as.matrix(svyby(~ Q8_mr, ~ REGION4, milconf_design, svymean, na.rm = TRUE))
dm_regn <- cbind(c(rep(5,4)), dm_regn)
dm_regn[,2] <- c("Northeast", "Midwest", "South", "West")

dm_incm <- as.matrix(svyby(~ Q8_mr, ~ income5, milconf_design, svymean, na.rm = TRUE))
dm_incm <- cbind(c(rep(2,5)), dm_incm)
dm_incm[,2] <- c("Q5", "Q4", "Q3", "Q2", "Q1")

dm_empl <- as.matrix(svyby(~ Q8_mr, ~ unemployed, milconf_design, svymean, na.rm = TRUE))
dm_empl <- cbind(c(rep(3,2)), dm_empl)
dm_empl[,2] <- c("Employed", "Unemployed")

dm_mart <- as.matrix(svyby(~ Q8_mr, ~ married, milconf_design, svymean, na.rm = TRUE))
dm_mart <- cbind(c(rep(1,2)), dm_mart)
dm_mart[,2] <- c("Single", "Married")

dm_geog <- as.matrix(svyby(~ Q8_mr, ~ URBAN3, milconf_design, svymean, na.rm = TRUE))
dm_geog <- cbind(c(rep(4,3)), dm_geog)
dm_geog[,2] <- c("Urban", "Suburban", "Rural")

dm_educ <- as.matrix(svyby(~ Q8_mr, ~ edulvl, milconf_design, svymean, na.rm = TRUE))
dm_educ <- cbind(c(rep(7,4)), dm_educ)
dm_educ[,2] <- c("HS or lower", "Some college", "College", "Graduate")

demo_preds1 <- as.data.frame(rbind(dm_mart, dm_incm, dm_empl,
                                   dm_geog, dm_regn, dm_genr,
                                   dm_educ, dm_race, dm_relg,
                                   dm_soco, dm_fami, dm_mvet,
                                   dm_gend, dm_ideo, dm_prty))
demo_preds1$link <- as.numeric(as.character(demo_preds1$Q8_mr)) * 100
demo_preds1$ci_lo <- (demo_preds1$link - 196 * as.numeric(as.character(demo_preds1$se))) 
demo_preds1$ci_hi <- (demo_preds1$link + 196 * as.numeric(as.character(demo_preds1$se)))
#demo_preds1$link2 <- as.factor(demo_preds1$link)

fig8_pts <- ggplot(demo_preds1, aes(y = link, x = V1, ymin = ci_lo, ymax = ci_hi, 
                                    label = married)) + 
  geom_linerange(colour = "gray65") +
  geom_point() + 
  geom_text_repel(size = 2.5, direction = "x", vjust = -1.3) +
  scale_x_discrete(breaks = c(1:15),
                   labels = c("Marital Status", "Income", "Employment",
                              "Geography", "Region", "Generation",
                              "Education", "Race", "Religion",
                              "Social Contact", "Family Contact", "Veteran Status",
                              "Gender", "Ideology", "Party ID")) +
  coord_flip() + xlab("") + ylab("Predicted Support for the Military (%)") + 
  ylim(40,90) + theme_bw()
ggsave("~/Dropbox/Public_Conf_Mil/Chapter 2 Figures/fig2-9_w1.jpeg", plot = fig8_pts, dpi = 700)

#fig8_op2 <- ggplot(demo_preds1, aes(y = link2, x = V1, label = party)) + 
#  geom_bar(position = position_dodge(width = 1, preserve = "single"), width = 0.35,
#           stat = "identity") + 
#  geom_text(size = 2, position = position_dodge(width = 1), hjust = -0.3) +
#  coord_flip() + xlab("") + ylab("Predicted Support for the Military") + 
#  scale_x_discrete(breaks = c(1:15),
#                   labels = c("Party ID", "Ideology", "Gender", "Race",
#                              "Religion", "Veteran", "Social Contact", "Family",
#                              "Generation", "Region", "Income", "Employment",
#                              "Martial Status", "Geography", "Education")) +
#  scale_y_discrete(breaks = c(40:90)) +
#  theme_bw()
#ggsave("demo_fig8-2.jpeg", plot = fig8_op2, dpi = 500)

#### FIGURE 2.9 - DEMOGRAPHIC PREDICTIONS ####

demo_mod <- svyglm(Q8_mr ~ dem + rep + ideo_m + male + black + hispanic + other_race + income5 +
                     edulvl + unemployed + boomer + genx + millen + activeduty + vet +
                     social + family + midwest + south + west + catholic + christian + norelig +
                     city + rural + married,
                   design = milconf_design)

### Predictions by Demographic Profiles

## Demo 1: Asian female student
demo1 <- data.frame(dem = 1, rep = 0, ideo_m = 2,
                    male = 0, black = 0, hispanic = 0, other_race = 1,
                    income5 = 2, edulvl = 3, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 0)
demo1pred <- as.data.frame(predict(demo_mod, newdata = demo1, se.fit = TRUE))

## Demo 2: Black low-income single female
demo2 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 0, black = 1, hispanic = 0, other_race = 0,
                    income5 = 2, edulvl = 2, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo2pred <- as.data.frame(predict(demo_mod, newdata = demo2, se.fit = TRUE))

## Demo 3: CA white collar male worker
demo3 <- data.frame(dem = 0, rep = 0, ideo_m = 3,
                    male = 1, black = 0, hispanic = 0, other_race = 1,
                    income5 = 5, edulvl = 4, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo3pred <- as.data.frame(predict(demo_mod, newdata = demo3, se.fit = TRUE))

## Demo 4: Hispanic male migrant worker
demo4 <- data.frame(dem = 0, rep = 0, ideo_m = 3,
                    male = 1, black = 0, hispanic = 1, other_race = 0,
                    income5 = 1, edulvl = 1, unemployed = 0,
                    boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 1)
demo4pred <- as.data.frame(predict(demo_mod, newdata = demo4, se.fit = TRUE))

## Demo 5: CA liberal male PhD
demo5 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 1, black = 0, hispanic = 0, other_race = 0,
                    income5 = 4, edulvl = 4, unemployed = 0,
                    boomer = 1, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo5pred <- as.data.frame(predict(demo_mod, newdata = demo5, se.fit = TRUE))

## Demo 6: NE Boomer woman
demo6 <- data.frame(dem = 1, rep = 0, ideo_m = 1,
                    male = 0, black = 0, hispanic = 0, other_race = 0,
                    income5 = 5, edulvl = 3, unemployed = 0,
                    boomer = 1, genx = 0, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 0,
                    midwest = 0, south = 0, west = 0,
                    catholic = 0, christian = 0, norelig = 1,
                    city = 0, rural = 0, married = 1)
demo6pred <- as.data.frame(predict(demo_mod, newdata = demo6, se.fit = TRUE))

## Demo 7: Midwestern soccer Mom
demo7 <- data.frame(dem = 0, rep = 0, ideo_m = 4,
                    male = 0, black = 0, hispanic = 0, other_race = 0,
                    income5 = 4, edulvl = 3, unemployed = 0,
                    boomer = 0, genx = 1, millen = 0,
                    activeduty = 0, vet = 0, social = 0, family = 1,
                    midwest = 1, south = 0, west = 0,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 0, married = 1)
demo7pred <- as.data.frame(predict(demo_mod, newdata = demo7, se.fit = TRUE))

## Demo 8: Black female NCO
demo8 <- data.frame(dem = 1, rep = 0, ideo_m = 2,
                    male = 0, black = 1, hispanic = 0, other_race = 0,
                    income5 = 3, edulvl = 2, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 1, vet = 1, social = 1, family = 1,
                    midwest = 0, south = 1, west = 0,
                    catholic = 0, christian = 1, norelig = 1,
                    city = 1, rural = 0, married = 0)
demo8pred <- as.data.frame(predict(demo_mod, newdata = demo8, se.fit = TRUE))

## Demo 9: Post-911 white male veteran
demo9 <- data.frame(dem = 0, rep = 1, ideo_m = 5,
                    male = 1, black = 0, hispanic = 0, other_race = 0,
                    income5 = 3, edulvl = 3, unemployed = 0,
                    boomer = 0, genx = 0, millen = 1,
                    activeduty = 0, vet = 1, social = 1, family = 0,
                    midwest = 0, south = 0, west = 1,
                    catholic = 1, christian = 0, norelig = 0,
                    city = 0, rural = 1, married = 1)
demo9pred <- as.data.frame(predict(demo_mod, newdata = demo9, se.fit = TRUE))

## Demo 10: Hispanic male marine
demo10 <- data.frame(dem = 0, rep = 0, ideo_m = 4,
                     male = 1, black = 0, hispanic = 1, other_race = 0,
                     income5 = 3, edulvl = 3, unemployed = 0,
                     boomer = 0, genx = 0, millen = 1,
                     activeduty = 1, vet = 1, social = 1, family = 1,
                     midwest = 0, south = 0, west = 1,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 1, rural = 0, married = 0)
demo10pred <- as.data.frame(predict(demo_mod, newdata = demo10, se.fit = TRUE))

## Demo 11: NASCAR Dad
demo11 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 3, edulvl = 3, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 1,
                     midwest = 0, south = 1, west = 0,
                     catholic = 0, christian = 0, norelig = 1,
                     city = 0, rural = 1, married = 0)
demo11pred <- as.data.frame(predict(demo_mod, newdata = demo11, se.fit = TRUE))

## Demo 12: NE white collar worker
demo12 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 4, edulvl = 4, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 0, social = 1, family = 0,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 1, rural = 0, married = 1)
demo12pred <- as.data.frame(predict(demo_mod, newdata = demo12, se.fit = TRUE))

## Demo 13: Male Midwestern farmer
demo13 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 2, edulvl = 2, unemployed = 0,
                     boomer = 0, genx = 1, millen = 0,
                     activeduty = 0, vet = 1, social = 1, family = 1,
                     midwest = 1, south = 0, west = 0,
                     catholic = 1, christian = 0, norelig = 0,
                     city = 0, rural = 1, married = 1)
demo13pred <- as.data.frame(predict(demo_mod, newdata = demo13, se.fit = TRUE))

## Demo 14: Female Fox News watcher
demo14 <- data.frame(dem = 0, rep = 1, ideo_m = 6,
                     male = 0, black = 0, hispanic = 0, other_race = 0,
                     income5 = 3, edulvl = 2, unemployed = 0,
                     boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 0, social = 0, family = 1,
                     midwest = 0, south = 0, west = 0,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 1)
demo14pred <- as.data.frame(predict(demo_mod, newdata = demo14, se.fit = TRUE))

## Demo 15: Male Republican Korean war veteran
demo15 <- data.frame(dem = 0, rep = 1, ideo_m = 7,
                     male = 1, black = 0, hispanic = 0, other_race = 0,
                     income5 = 4, edulvl = 3, unemployed = 0,
                     boomer = 1, genx = 0, millen = 0,
                     activeduty = 0, vet = 1, social = 0, family = 0,
                     midwest = 0, south = 0, west = 1,
                     catholic = 0, christian = 1, norelig = 0,
                     city = 0, rural = 0, married = 0)
demo15pred <- as.data.frame(predict(demo_mod, newdata = demo15, se.fit = TRUE))

demo_preds <- rbind(demo1pred, demo2pred, demo3pred, demo4pred, demo5pred,
                    demo6pred, demo7pred, demo9pred, demo8pred, demo11pred, 
                    demo10pred, demo12pred, demo14pred, demo13pred, demo15pred)
demo_preds$ci_lo <- demo_preds$link - 1.96*demo_preds$SE
demo_preds$ci_hi <- demo_preds$link + 1.96*demo_preds$SE
demo_preds$count <- c(1:15)

demo_plot <- ggplot(demo_preds) +
  geom_pointrange(aes(x = count, y = link, ymin = ci_lo, ymax = ci_hi)) +
  coord_flip() + xlab("") + ylab("Predicted Support for the Military") +
  scale_x_continuous(breaks = c(1:15),
                     labels = c("Asian female student", "Black low-income single female",
                                "CA male white collar worker", 
                                "Hispanic male migrant worker",
                                "CA liberal male PhD", "NE Boomer woman", 
                                "Midwestern soccer mom", "Post-9/11 white male veteran",
                                "Black female NCO", 
                                "NASCAR Dad", "Hispanic male marine",
                                "NE white collar worker", 
                                "Female Fox News watcher",
                                "Male Midwestern farmer",
                                "Male Republican Korean war veteran")) +
  theme_bw()
demo_plot
ggsave("fig2-10_w1.jpeg", plot = demo_plot, dpi = 1000)
